{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Modeling transmon qubit Cooper-pair box Hamiltonian in the charge basis "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "(Zlatko Minev, Christopher Warren, Nick Lanzillo 2021)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This module models the transmon qubit in the cooper-pair charge basis, assuming wrapped junction phase variable. The Hamiltonian is given by: \n",
    "\n",
    "$$\n",
    "\\hat{H}=4E_{C}\\left(\\hat{n}-n_{g}\\right)-E_{J}\\cos\\left(\\hat{\\phi}\\right)\\,,\n",
    "$$\n",
    "\n",
    "where $E_{C}$ is the charging energy, $E_{J}$ is the Josephson energy, $\\hat n$ is the number of Cooper pairs transferred between charge islands, $\\hat{\\phi}$ is the gauge-invariant phase difference between charge islands, and $n_{g}$ is effective offset charge of the device. Expressions for the charging energy, Josephson energy and offset charge can be written as:\n",
    "\n",
    "$$\n",
    "E_{C}=\\frac{e^{2}}{2C_{\\Sigma}}\\,,\\qquad n_{g}=-\\frac{C_{d}\\dot{\\Phi}_{s}\\left(t\\right)}{2e}\\:,\\qquad E_{J}=\\frac{\\phi_{0}^{2}}{L_{J}}\\,,\n",
    "$$\n",
    "\n",
    "where $C_{\\Sigma} = C_{J}+C_{B}+C_{g}$ (the sum of the Josephson capacitance, shunting capacitance and gate capacitance), $L_{J}$ is the inductance of the Josephson junction, and $\\phi$ is the magnetic flux. \n",
    "\n",
    "The variables are\n",
    "$$\n",
    "\\hat{\\phi}\\equiv\\frac{\\hat{\\Phi}}{\\phi_{0}},\\qquad\\hat{n}\\equiv\\frac{\\hat{Q}}{2e}\\,,\n",
    "$$\n",
    "\n",
    "Observe that $\\hat \\phi$ and $\\hat n$ are both dimensiuonless, and they obey the commutation relationship:\n",
    "\n",
    "$$\n",
    "[\\hat{\\phi}, \\hat{n}] = i\n",
    "$$\n",
    "\n",
    "\n",
    "The Hamiltonian can be written in the charge ($\\hat n$) basis as: \n",
    "\n",
    "$$H=4E_\\text{C}(\\hat{n}-n_g)^2-\\frac{1}{2}E_\\text{J}\\sum_n(|n\\rangle\\langle n+1|+\\text{h.c.}),$$\n",
    "Where $\\hat{n} = \\sum_{n=-\\inf}^{\\inf} |n\\rangle\\langle n|$\n",
    "\n",
    "### Hcpb class\n",
    "\n",
    "Hamiltonian-model Cooper pair box (Hcpb) class.\n",
    "\n",
    "Used to model analytically the CPB Hamiltonian quickly\n",
    "and efficiently. Solves in charge basis tridiagonal eigenvalue\n",
    "problem for arbitrary Ej, Ec, ng values.\n",
    "\n",
    "As long as nlevels remains fixed the number of charge states\n",
    "considered does not change and it does not recreate the arrays,\n",
    "just recomputes the properties\n",
    "\n",
    "Returns all properties of interest for the CPB.\n",
    "\n",
    "This model is closer to the analytic solution than the Duffing oscillator model.\n",
    "Can work backwards from target qubit parameters to get the Ej, Ec or use\n",
    "input Ej, Ec to find the spectrum of the Cooper Pair Box.\n",
    "\n",
    "    @author: Christopher Warren (Chalmers University of Technology), updated by Zlatko K. Minev (IBM Quantum)\n",
    "    @date: 2020, 2021\n",
    "\n",
    "\n"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's start by importing the key modules for this demo: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "%load_ext autoreload\r\n",
    "%autoreload 2\r\n",
    "%config IPCompleter.greedy = True\r\n",
    "%matplotlib inline\r\n",
    "%config Completer.use_jedi = False\r\n",
    "%config InlineBackend.figure_format = 'svg'\r\n",
    "\r\n",
    "from qiskit_metal.analyses.hamiltonian.transmon_charge_basis import Hcpb\r\n",
    "from qiskit_metal.analyses.hamiltonian.transmon_CPB_analytic import Hcpb_analytic\r\n",
    "import pandas as pd\r\n",
    "import matplotlib.pyplot as plt\r\n",
    "import numpy as np"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Let's model a transmon "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Energy levels"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can get a feel for how to use the Hcpb class by plotting the first few eigenvalues as a function of offset charge, similar to the plots reported in Phys. Rev. A 76, 042319 (2007.) Let's start by defining the range of offset charge from -2.0 to +2.0 and also by defining a normalization for the eigenvalues, which will be the transition energy between the first two states evaluated at the degenercy point where ng=0.5. For this exercise, we'll take Josephson Energy to be equal to the charging energy:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x = np.linspace(-2.0,2.0,101) # this represents the charging energy (ng)\r\n",
    "H_norm = Hcpb(nlevels=2, Ej=1000.0, Ec=1000.0, ng=0.5) # Hamiltonian definition \r\n",
    "norm = H_norm.fij(0,1) # normalization constant "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we'll sweep over the offset charge and calculate the first three eigenvalues for a given value of ng. We'll need to define a new Hamiltonian for this. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "E0, E1, E2 = [], [], []\r\n",
    "\r\n",
    "# For a given value of offset charge (ng, represented by x) we will calculate the CPB Hamiltonian using the previously assigned values of E_J and E_C. Then we calculate the eigenvalue for a given value of m.\r\n",
    "for i in x: \r\n",
    "    H = Hcpb(nlevels=3, Ej=1000.0, Ec=1000.0, ng=i)\r\n",
    "    E0.append(H.evalue_k(0)/norm)\r\n",
    "    E1.append(H.evalue_k(1)/norm)\r\n",
    "    E2.append(H.evalue_k(2)/norm)\r\n",
    "\r\n",
    "# define the minimum of E0 and set this to E=0\r\n",
    "floor = min(E0) \r\n",
    " \r\n",
    "plt.plot(x, E0 - floor, 'k', label=\"m=0\")\r\n",
    "plt.plot(x, E1 - floor, 'r', label=\"m=1\")\r\n",
    "plt.plot(x, E2 - floor, 'b', label=\"m=2\")\r\n",
    "plt.xlabel(\"ng\")\r\n",
    "plt.ylabel(\"Em/E01\")\r\n",
    "plt.legend(title=\"Energy Level:\", loc='upper right')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Comparing with the Analytic Expressions for Energy"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can compare calculated eigenvalues with the analytic solutions by using the \"Hcpb_analytic\" class, which calculates the transmon eigenvalues analytically using Mathieu characteristic values instead of a matrix-based approach. Let's compare the calculated values of the lowest energy at zero offset charge in both cases:  "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# this is using the Hcpb approach as above, solving the charge basis tridiagonal eigenvalue problem:\r\n",
    "H_CPB = Hcpb(nlevels=15, Ej=13971.3, Ec=295.2, ng=0.0)\r\n",
    "\r\n",
    "# this using the Hcpb_analytic class, which solves using the exact (analytic) solutions in terms of Mathieu characteristic values: \r\n",
    "H_CPB_analytic = Hcpb_analytic(Ej=13971.3, Ec=295.2, ng=0.0)\r\n",
    "\r\n",
    "# print and compare energies \r\n",
    "print(\"E0 (HCPB):\", H_CPB.evalue_k(0))\r\n",
    "print(\"E0 (HCPB analytic):\", H_CPB_analytic.evalue_k(0))\r\n",
    "print(\"Error:\", 100*(H_CPB_analytic.evalue_k(0) - H_CPB.evalue_k(0)) / H_CPB_analytic.evalue_k(0))\r\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As we can see above, the calculated eigenvalues using the Hcbp class match the analytic values extremely well. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Wavefunctions "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's define a new Hamiltonian, this time with $E_{J} >> E_{C}$ and an offset charge of $ng=0.001$. We can calculate the transition energy between the lowest two states as well as the anharmonicity with the following: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "H = Hcpb(nlevels=3, Ej=13971.3, Ec=295.2, ng=0.001)\r\n",
    "print(f\"\"\"\r\n",
    "Transmon frequencies \r\n",
    "\r\n",
    " ω01/2π = {H.fij(0,1): 6.0f} MHz\r\n",
    "   α/2π = {H.anharm(): 6.0f} MHz\r\n",
    "\"\"\")\r\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that both the transition energy and the anharmonicity are read out in units of Megahertz (MHz.) "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can plot the eigenstates (wavefunctions) of the transmon qubit using the commands below: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import matplotlib.pyplot as plt\r\n",
    "for k in range (3):\r\n",
    "    ψ, θ = H.psi_k(k,101)\r\n",
    "    plt.plot(θ, ψ.real+ψ.imag, label=f\"|{k}>\") # it's in either quadrature, but not both\r\n",
    "plt.xlabel(\"Junction phase θ (wrapped in the interval [-π, π])\")\r\n",
    "plt.ylabel(\"Re(ψ(θ))\")\r\n",
    "plt.legend(title=\"Level\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Verifying Orthonormality of the Wavefunctions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we can verify the orthonormality of the wavefunctions. Let's first take the first two eigenstates and verify that their inner product is zero, thereby confirming orthogonality. Note that since the wavefunctions can be complex, we need to take the complex conjugate of $\\Psi1$. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "Psi0, theta0 = H.psi_k(0)\r\n",
    "Psi1, theta1 = H.psi_k(1)\r\n",
    "print(np.dot(Psi0,Psi1.conj()))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Indeed, we see that the dot product is essentially zero (within numerical precision.) Next, let's take the inner product of the first eigenstate with itself, checking that we get an output of unity:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "print(np.dot(Psi0, Psi0.conj()))\r\n",
    "print(np.dot(Psi1, Psi1.conj()))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Indeed we see that the dot products are essentially equal to unity, confirming that the states are appropriately normalized.  "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Additional Analysis: Charge Dispersion, Energy Level Differences, Anharmonicity and Dephasing Time (T2)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Charge Dispersion"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The peak-to-peak value of the charge dispersion for the mth energy level is given by the expression: $\\epsilon_{m} = E_{m}(n_{g}=0.5) - E_{m}(n_{g}=0.0)$. We can plot $\\epsilon_{m}/E_{01}$ as a function of $E_{J}/E_{C}$ for the first few energy levels and reproduce the figure published in Phys. Rev. A 76, 042319 (2007) (Figure 4(a)). \n",
    "\n",
    "We can start by defining a value of charging energy and creating empty lists for $\\epsilon_{0}$ through $\\epsilon_{4}$: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "E_c=100.0 # charging energy \r\n",
    "epsilon0, epsilon1, epsilon2, epsilon3 = [], [], [], []    # charge dispersion for m=0 through m=4\r\n",
    "x = np.linspace(1,140,101)           # this this ratio of Ej/Ec which will go on the x-axis. "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we simply evaluate the expression given above for $\\epsilon_{m}$ based on $E_{m}$ and $E_{0}$. We use two separate Hamiltonians to do this; one evaluated at $n_{g}=0.5$ and one evaluated at $n_{g}=0.0$. We also normalize by the transition energy between the lowest two states evaluated at the degeneracy point ($E_{01}$.) Finally, we populate the lists each $\\epsilon_{m}$. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "for i in x:\r\n",
    "    E_j = i*E_c \r\n",
    "    H_zero = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.0)\r\n",
    "    H_half = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    \r\n",
    "    H_norm = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    norm = H_norm.fij(0,1)                         # normalization constant \r\n",
    "    \r\n",
    "    epsilon0.append(abs(H_half.evalue_k(0) - H_zero.evalue_k(0))/norm)\r\n",
    "    epsilon1.append(abs(H_half.evalue_k(1) - H_zero.evalue_k(1))/norm)\r\n",
    "    epsilon2.append(abs(H_half.evalue_k(2) - H_zero.evalue_k(2))/norm)\r\n",
    "    epsilon3.append(abs(H_half.evalue_k(3) - H_zero.evalue_k(3))/norm)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can plot these values to see the exponential decrease in the charge dispersion with increasing $E_{J}/E_{C}$: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.plot(x, epsilon0, 'k', label=\"m=0\")\r\n",
    "plt.plot(x, epsilon1, 'b', label=\"m=1\")\r\n",
    "plt.plot(x, epsilon2, 'r', label=\"m=2\")\r\n",
    "plt.plot(x, epsilon3, 'g', label=\"m=3\") \r\n",
    "plt.yscale(\"log\")\r\n",
    "plt.xlabel(\"EJ/Ec\")\r\n",
    "plt.ylabel(\"Epsilon/E01\")\r\n",
    "plt.legend(title=\"Energy Level\", loc='upper right')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Energy Level Differences"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also evaluate the energy level difference ($E_{m0} = E_{m} - E_{0}$) evaluated at the degeneracy point ($n_{g}=0.5$) and compare to Fig. 4(b) of Phys. Rev. A 76, 042319 (2007). To do this, we just need to create empty lists for the energy levels ($E_{0}$ through $E_{3}$) as well as the energy level differences ($E_{00}$ through $E_{30}$.) "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "E00, E10, E20, E30 = [], [], [], [] \r\n",
    "E0, E1, E2, E3 = [], [], [], [] "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we'll sweep over $E_{J}/E_{C}$ from 0 to 140 (using the variable x defined above) and at each point, we'll construct the Hamiltonian and take the difference in eigenvalues evaluated at the degeneracy point $n_{g}=0.5$. We'll also normalize the energy level differences by the charging energy to be consistent with Fig. 4(b) in the above reference: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "for i in x:\r\n",
    "    H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    E0 = H.evalue_k(0) \r\n",
    "    E1 = H.evalue_k(1)\r\n",
    "    E2 = H.evalue_k(2)\r\n",
    "    E3 = H.evalue_k(3) \r\n",
    "    E00.append((E0 - E0)/E_c)\r\n",
    "    E10.append((E1 - E0)/E_c)\r\n",
    "    E20.append((E2 - E0)/E_c)\r\n",
    "    E30.append((E3 - E0)/E_c)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can plot these results and see how the energy level differences increase with increasing $E_{J}/E_{C}$ ratio: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.plot(x,E00,'k',label=\"m0\")\r\n",
    "plt.plot(x,E10,'b',label=\"m=1\")\r\n",
    "plt.plot(x,E20,'g',label=\"m=2\")\r\n",
    "plt.plot(x,E30,'r',label=\"m=3\") \r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"Em0/Ec\")\r\n",
    "plt.legend(title=\"Energy Levels\", loc='upper right')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Anharmonicity"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We know that for the transmon qubit, having the Josephson Energy much larger than the charging energy ($E_{J} >> E_{C}$) results in a decrease in anharmonicity. The latter is critical for a functional qubit in which the energy difference between the lowest two states ($E_{01}$) is sufficiently different than the energy difference between the second and third states, $E_{12}$. The absolute anharmonicity is defined as $\\alpha = E_{12} - E_{01}$, while the relative anharmonicity is defined as $\\alpha_{r} = \\alpha/E_{01}$. \n",
    "\n",
    "We can easily make a plot of the anharmonicity as a function of $E_{J}/E_{C}$ using the Hcpb class. Let's have the ratio of $E_{J}/E_{C}$ (which we'll call x) vary from 0 to 80, and then we'll create empty lists for $\\alpha$ and $\\alpha_r$. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x = np.linspace(0,80,101)   #EJ/EC\r\n",
    "alpha = [] \r\n",
    "alpha_r = [] "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we'll just sweep over x and at each value we'll calculate both the absolute and relative anharmonicity. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "for i in x:\r\n",
    "    H_anharm = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    alpha.append(H_anharm.anharm()) \r\n",
    "    alpha_r.append(H_anharm.anharm()/H_anharm.fij(0,1))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note we get a warning here because the relative anharmonicity blows up as Ej/Ec goes to zero. Then we can plot the results:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.figure(1)\r\n",
    "plt.subplot(131)\r\n",
    "plt.plot(x,alpha)\r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"alpha\") \r\n",
    "plt.subplot(133)\r\n",
    "plt.plot(x,alpha_r)\r\n",
    "plt.ylim(-0.2, 1.0) \r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"alpha_r\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Indeed we see that the anharmonicity decays with the inverse of $E_{J}/E_{C}$ for small values of $E_{J}/E_{C}$ before reaching a minimum just before $E_{J}/E_{C} \\approx 20.0$, then changing sign and approaching zero as $E_{J}/E_{C}$ approaches infinity. This matches very closely to the results found in Figure 5 of Phys. Rev. A 76, 042319 (2007). "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Dephasing Time (T2)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can estimate the qubit dephasing time (T2) due to charge noise by the following expression: $T_{2} = \\frac{\\hbar}{A\\pi |\\epsilon_{1}|}$ where $A$ is on the order of $1E-4$ according to Phys. Rev. A 76, 042319 (2007). Since this is essentially just the inverse of the charge dispersion for $\\epsilon_{1}$, we can easily calculate T2 as a function of $E_{J}/E_{C}$ with the following:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x = np.linspace(0,80,101)     # ratio of Ej/Ec, varying from 0 to 80\r\n",
    "T2 = []                       # empty list for T2 "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we'll just calculate the T2 time as a function of $E_{J}/E_C{}$, noting that the output of the Hcpb eigenvalue caluclation is in units of E/h ~ MHz (1E6 Hz). So our calculation simplifies to: $T_{2} = \\frac{1.0}{2*(1E-4)*(1E6) |\\epsilon_{1}|}$ where $\\epsilon_{1}$ is the direct output of the Hcpb eigenvalue calculation: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "E_c = 1000.0 \r\n",
    "for i in x: \r\n",
    "    H_half = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    H_zero = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.0)\r\n",
    "    eps = abs(H_half.evalue_k(1) - H_zero.evalue_k(1))    \r\n",
    "    T2.append(1.0/(2.0*(1E-4)*(1E6)*eps) )\r\n",
    "\r\n",
    "plt.plot(x, T2)\r\n",
    "plt.yscale(\"log\")\r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"T2 (sec)\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Indeed, we see that the dephasing time increases exponentially with increasing $E_{J}/E_{C}$, which is one of the critical features of the transmon qubit. By increasing the $E_{J}/E_{C}$ ratio, we reduce sensitivity to charge noise without sacrificing too much anharmonicity, resulting in greatly improved dephasing time. This plot matches very closely to Fig. 5(c) in Phys. Rev. A 76, 042319 (2007)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Qutip simulation "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Diagonal Hamiltonian\n",
    "\n",
    "Wrapper around Qutip to output the diagonalized\n",
    "Hamiltonian truncated up to n levels of the transmon\n",
    "for modeling"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "H.h0_to_qutip(6)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Coupling and number operator\n",
    "Wrapper around Qutip to output the number operator (charge)\n",
    "for the Transmon Hamiltonian in the energy eigen-basis.\n",
    "Used for computing the\n",
    "coupling between other elements in the system."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "H.n_to_qutip(6)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can use these matrices to perform some additional analysis. As an example, we can plot the fluctuation in the number of Cooper pairs as a function of $E_{J} / E_{C}$ for the first three energy levels. This requires calculating the variance in the number of Cooper pairs, defined as: $\\sqrt{<n^2> - <n>^2 }$ where $n$ is calculated for a given energy level m. \n",
    "\n",
    "So let's first set our x-axis (representing $E_{J}/E_{C}$) to run from 1 to 120 and we'll create empty lists to represent the variances for the m=0, m=1, and m=2 states. We'll also define a value of $E_{C}$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x = np.linspace(1,120, 101)\r\n",
    "E_c=350.0\r\n",
    "var0 = [] \r\n",
    "var1 = []\r\n",
    "var2 = [] "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we'll define the Hamiltonian matrix and calculate both $<n>$ and $<n^2>$. Then we'll take the square root of the difference between the two for each value of $E_{J}/E_{C}$. Lastly, these values get appended to the empty lists that we defined above: "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "for i in x:\r\n",
    "    H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    n_squared = H.n_to_qutip(6)*H.n_to_qutip(6)\r\n",
    "    n = H.n_to_qutip(6)\r\n",
    "    fluc0 = np.sqrt(np.real(n_squared[0,0] - (n[0,0])*(n[0,0])))\r\n",
    "    fluc1 = np.sqrt(np.real(n_squared[1,1] - (n[1,1])*(n[1,1])))\r\n",
    "    fluc2 = np.sqrt(np.real(n_squared[2,2] - (n[2,2])*(n[2,2])))\r\n",
    "    var0.append(fluc0)\r\n",
    "    var1.append(fluc1)\r\n",
    "    var2.append(fluc2)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can plot these and check that the resulting curves match quite closely to those found in Fig. 6 of Koch et al. (2007) PRA:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.plot(x,var0,'k', label=\"m=0\")\r\n",
    "plt.plot(x,var1,'r', label=\"m=1\")\r\n",
    "plt.plot(x,var2,'b', label=\"m=2\")\r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"Variance\")\r\n",
    "plt.legend(title=\"Energy Levels\", loc='upper right')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Physically, the above plot demonstrates that the variations in the number of Cooper pairs is on the order of 1-2 for gorund and first excited states of the CPB operated in the transmon regime."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also investigate how the off-diagonal matrix elements of the CPB Hamiltonian evolve as a function of the $E_{J}/E_{C}$ ratio. Let's create four empty lists ($m_{10}$, $m_{20}$, $m_{30}$ and $m_{21}$) where $m_{ij} = <i | n | j>$. Then we can loop over the ratio of $E_{J}/E_{C}$ and at each value, we calculate the off-diagonal matrix element:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "m10 = [] \r\n",
    "m20 = [] \r\n",
    "m30 = [] \r\n",
    "m21 = [] \r\n",
    "\r\n",
    "# For a given value of offset charge (ng, represented by x) we will calculate the CPB Hamiltonian using the previously assigned values of E_J and E_C. Then we calculate the eigenvalue for a given value of m.\r\n",
    "for i in x: \r\n",
    "    H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)\r\n",
    "    m10.append(np.real(H.n_to_qutip(6)[1,0]))\r\n",
    "    m20.append(np.real(H.n_to_qutip(6)[2,0]))\r\n",
    "    m30.append(np.real(H.n_to_qutip(6)[3,0]))\r\n",
    "    m21.append(np.real(H.n_to_qutip(6)[2,1])) "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plotting these off-diagonal matrix elements is then straightforward:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.plot(x,m10,'k',label=\"<1|n|0>\")\r\n",
    "plt.plot(x,m20,'r',label=\"<2|n|0>\")\r\n",
    "plt.plot(x,m30,'y',label=\"<3|n|0>\")\r\n",
    "plt.plot(x,m21,'b',label=\"<2|n|1>\")\r\n",
    "plt.xlabel(\"Ej/Ec\")\r\n",
    "plt.ylabel(\"<i|n|j>\")\r\n",
    "plt.legend(title=\"Energy Levels\", loc='upper right')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We see that the above plot closely matches Fig. 7 of Koch et al. (2007) PRA. Physically, we see that the coupling between adjacent transmon states ($<1|n|0>$ and $<2|n|1>$) is the only significant coupling in the $E_{J} >> E_{C}$ regime. The coupling between non-adjacent states (like m=0 and m=2) is essentially zero at all values of $E_{J}/E_{C}$. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Experimental  "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's use the \"params_from_spectrum\" function to calculate the target Ej and Ec values for a desired transmon frequency and anharmonicity. We'll use values of Ej=13971.3 MHz and Ec=295.2 MHz. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# 13971.3, Ec=295.2\r\n",
    "ω, α = 5431, -341\r\n",
    "EjEc = H.params_from_spectrum(ω, α) # set self.Ej, Cj\r\n",
    "print(EjEc)\r\n",
    "print(\"transmon frequency:\", H.fij(0,1), \"anharmonicity:\", H.anharm())\r\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also calculate the value of Ej given the value of Ec and the transmon frequency:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "Ej = H.params_from_freq_fixEC(ω, 295.17)\r\n",
    "print(\"Ej:\", Ej)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We see that this is close to the starting vaue of 13971.3 MHz that we began with. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# New section on integrating sc_qubits. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This section is in development. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import scqubits as scq"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "qubit = scq.Transmon(\r\n",
    "    EJ=13.97124102543,\r\n",
    "    EC=0.295179,\r\n",
    "    ng=0.001,\r\n",
    "    ncut=40,\r\n",
    "    truncated_dim=4     # after diagonalization, we will keep 3 levels\r\n",
    ")\r\n",
    "evals = qubit.eigenvals(evals_count=12)\r\n",
    "print(f\"freq = {(evals[1] - evals[0])* 1000:.0f} MHz\")\r\n",
    "print(f\"alpha = {((evals[2] - evals[1]) - (evals[1] - evals[0]))* 1000:.0f} MHz\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "qubit.plot_n_wavefunction()\r\n",
    "qubit.plot_phi_wavefunction(which=[0,1,2], mode='real')\r\n",
    "qubit.hamiltonian()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import numpy as np\r\n",
    "ng_list = np.linspace(-2, 2, 220)\r\n",
    "qubit.plot_evals_vs_paramvals('ng', ng_list, evals_count=6, subtract_ground=False);"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Transmon and the Oscillator"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This section is currently in development. "
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}